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THEORY  OF  THE  LATTICE  BOLTZMANN  METHOD:  LATTICE  BOLTZMANN 

MODELS  FOR  NON-IDEAL  GASES 

LI-SHI  LUO* 

Abstract.  In  this  paper  a  procedure  for  systematic  a  priori  derivation  of  the  lattice  Boltzmann  models  for 
non-ideal  gases  from  the  Enskog  equation  (the  modified  Boltzmann  equation  for  dense  gases)  is  presented. 
This  treatment  provides  a  unified  theory  of  lattice  Boltzmann  models  for  non-ideal  gases.  The  lattice 
Boltzmann  equation  is  systematically  obtained  by  discretizing  the  Enskog  equation  in  phase  space  and 
time.  The  lattice  Boltzmann  model  derived  in  this  paper  is  thermodynamically  consistent  up  to  the  order 
of  discretization  error.  Existing  lattice  Boltzmann  models  for  non-ideal  gases  are  analyzed  and  compared 
in  detail.  Evaluation  of  these  models  are  made  in  light  of  the  general  procedure  to  construct  the  lattice 
Boltzmann  model  for  non-ideal  gases  presented  in  this  work. 

Key  words,  kinetic  method,  Enskog  theory,  non-ideal  gases,  lattice  Boltzmann  equation 

Subject  classification.  Fluid  Mechanics 

1.  Introduction.  In  recent  years,  there  has  been  significant  progress  made  in  the  development  of  the 
lattice  Boltzmann  equation  (LBE)  method  [48,  29,  28,  6,  10,  1,  38],  a  novel  technique  developed  for  modeling 
various  complex  systems,  especially  complex  fluids.  One  particular  application  of  the  lattice  Boltzmann 
method  which  has  attracted  considerable  attention  is  the  modeling  of  inhomogeneous  fluids,  such  as  non¬ 
ideal  gases  or  multi-component  fluids  [21,  18,  19,  20,  59,  59,  61,  62,  51,  50,  27,  46].  These  flows  are  important, 
but  are  difficult  to  simulate  by  conventional  techniques  for  solving  the  Navier-Stokes  equations.  The  main 
difficulty  conventional  techniques  face  is  due  to  the  interfaces  in  inhomogeneous  flow.  Computationally, 
one  might  be  able  to  track  a  few,  but  certainly  not  very  many  interfaces  in  a  system.  It  is  therefore 
impractical  to  simulate  a  realistic  system,  which  is  inhomogeneous  in  density  or  composition,  by  directly 
solving  the  Navier-Stokes  equations  without  making  some  drastic  approximations.  One  can  also  view  this 
problem  from  a  different  perspective:  Interfaces  between  different  components  or  phases  of  a  fluid  system  are 
thermodynamic  effects  resulting  from  interactions  among  molecules.  To  solve  the  Navier-Stokes  equations, 
one  needs  to  know  the  equation  of  state,  which  is  usually  unknown  in  the  interfacial  regions.  It  is  therefore 
difficult  to  incorporate  thermodynamics  into  the  Navier-Stokes  equations  in  a  consistent  fashion,  although  the 
interfaces  are  precisely  the  result  due  to  thermodynamical  effects.  Hence  one  encounters  some  fundamental 
difficulties. 

There  exists  ample  evidence  that  models  based  on  the  lattice  Boltzmann  equation,  and  its  predecessor, 
the  lattice-gas  automata  (LG A)  [14,  65,  15],  and  other  gas  kinetic  models  [53,  68,  37],  are  particularly  suitable 
for  complex  systems  such  as  non-ideal  gases  and  multi-component  fluids  [59,  60,  61,  62,  51,  50,  27].  There  may 
be  profound  reasons  for  the  success  of  the  LGA  and  LBE  models  in  simulating  those  complex  systems.  The 
LGA  and  LBE  models  do  not  start  at  the  macroscopic  level;  instead,  they  start  at  mesoscopic  level  at  which 
one  can  use  a  “ potential ”  to  model  interactions  in  the  system.  Macroscopic  or  hydrodynamic  behaviors  of  the 
system  naturally  emerge  from  mesoscopic  dynamics,  provided  that  the  mesoscopic  dynamics  possesses  the 
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(email  address:  luo@icase.edu).  This  research  was  supported  by  the  National  Aeronautics  and  Space  Administration  under 
NASA  Contract  No.  NAS1-97046  while  the  author  was  in  residence  at  ICASE,  NASA  Langley  Research  Center,  Hampton,  VA 
23681-2199. 
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necessary  and  correct  conservation  laws  with  associated  symmetries  such  as  rotational  invariance,  Galilean 
invariance,  etc.  It  is  well  known  that  the  macroscopic  behavior  of  a  hydrodynamic  system  is  rather  insensitive 
to  the  microscopic  or  mesoscopic  details  —  the  details  of  microscopic  or  mesoscopic  dynamics  only  affect  the 
numerical  values  of  the  transport  coefficients.  This  observation  is  a  key  physical  insight  in  the  construction 
of  simplistic  kinetic  models  such  as  the  lattice  gas  automata  and  the  lattice  Boltzmann  equation. 

Historically,  the  lattice  Boltzmann  equation  was  first  developed  empirically  [48,  29,  28,  6]  from  its 
predecessor  —  the  lattice-gas  automata  [14,  65,  15].  This  empiricism  influences  even  the  most  recent  lattice 
Boltzmann  models  [59,  60,  62,  61,  51,  50,  27].  Empirical  lattice  Boltzmann  models  usually  have  some  inherent 
artifacts  which  are  not  yet  fully  understood.  One  particular  problem  with  non-ideal  gases  or  multi-component 
lattice  Boltzmann  models  is  the  thermodynamic  inconsistency:  The  so-called  “equilibrium  state”  in  these 
models  cannot  be  described  by  thermodynamics.  In  particular,  one  has  difficulties  in  defining  an  entropy  of 
the  system  systematically,  and  thus  leading  to,  for  instance,  the  inconsistency  between  the  thermodynamic 
pressure  and  the  kinetic  one  [9].  Although  this  issue  has  been  previously  mentioned  [62,  61],  no  progress  has 
been  made  in  solving  this  problem,  despite  its  paramount  importance. 

It  is  well  understood  that  the  original  Boltzmann  equation  only  describes  rarefied  gases;  it  does  not 
describe  dense  gases  or  liquids.  In  the  Boltzmann  gas  limit  (BGL),  N  — >•  oo,  m  — >•  0,  and  ro  — >•  0,  Nm  — >• 
finite,  Nr$  — >  finite,  and  Nr$  — >  0,  where  TV,  m,  and  ro  are  the  particle  number,  particle  mass,  and 
interaction  range,  respectively.  Thus,  in  the  BGL,  the  mean  free  path  l  ~  1/Nr^  remains  constant,  while 
the  total  interaction  volume  Nr%  goes  to  zero.  Therefore,  in  the  strict  thermodynamic  sense,  the  Boltzmann 
equation  only  retains  the  thermodynamic  properties  of  a  perfect  gas  —  there  is  no  contribution  to  the 
transport  of  molecular  properties  from  inter-particle  forces,  although  collisions  influenced  by  interparticle 
interaction  are  considered.  In  order  to  properly  describe  non-ideal  dense  gases,  the  effect  of  finite  particle 
size,  for  instance,  must  be  explicitly  considered.  It  was  Enskog  who  first  extended  the  Boltzmann  equation 
to  dense  gases  by  including  the  volume  exclusion  effect  along  the  rationale  of  van  der  Waals  theory  [5] ,  which 
leads  to  a  non-ideal  gas  equation  of  state.  The  Enskog  equation  (the  modified  Boltzmann  equation  for  dense 
gases)  can  indeed  describe  dense  gases  or  liquids  with  non-ideal  gas  equation  of  state  to  certain  extend.  The 
Enskog  equation  describes  a  system  consisting  of  hard  spheres  and  it  has  been  shown  that  the  hard-sphere 
system  captures  most  qualitative  properties  of  a  simple  liquid  [22,  57].  Furthermore,  the  revised  Enskog 
theory  seems  to  be  valid  for  a  wide  range  of  density  covering  gas,  liquid,  and  even  solid  [12]. 

It  has  been  recently  demonstrated  [25,  26]  that  the  lattice  Boltzmann  equation  can  be  directly  derived 
from  the  continuous  Boltzmann  equation.  The  method  proposed  in  Refs.  [25,  26]  is  a  general  procedure  to 
construct  the  lattice  Boltzmann  models  in  a  systematic  and  a  priori  fashion.  Through  this  procedure  we  can 
better  understand  the  approximation  made  in  the  lattice  Boltzmann  equation.  The  method  also  provides 
a  means  to  analyze  the  existing  lattice  Boltzmann  models.  In  this  paper,  the  method  of  Refs.  [25,  26]  is 
applied  to  obtain  the  lattice  Boltzmann  equation  for  non-ideal  gases  (which  have  non-ideal  gas  equation  of 
state).  The  lattice  Boltzmann  equation  for  non-ideal  gases  is  derived  from  the  Enskog  equation  for  dense 
gases.  The  obtained  lattice  Boltzmann  model  for  isothermal  non-ideal  gases  has  thermodynamic  consistency 
in  the  sense  of  approximation,  i.e.,  it  is  only  correct  up  to  the  order  of  discretization.  In  comparing  all 
existing  models,  we  would  like  to  stress  the  fact  that  many  defects  of  the  existing  LBE  models  are  due  to 
errors  made  at  the  level  of  fundamental  concepts,  rather  than  at  the  level  of  numerical  implementation. 

This  paper  is  a  detailed  extension  of  a  work  previously  published  [43],  and  is  organize  as  follows.  In  Sec.  2 
the  Enskog  equation  for  dense  gases  with  BGK  approximation  is  briefly  discussed.  In  Sec.  3  the  discretization 
procedure  to  obtain  the  lattice  Boltzmann  equation  for  non-ideal  gases  from  the  Enskog  equation  is  described. 
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The  discretization  in  time  and  phase  space,  small  velocity  expansion  of  the  equilibrium  distribution,  and 
realization  of  the  forcing  term  in  the  lattice  Boltzmann  equation  are  also  discussed  in  detail.  In  Sec.  4  the 
hydrodynamics  and  some  related  thermodynamical  quantities  of  the  model  are  given.  In  Sec.  5  the  model 
derived  in  this  work  is  compared  with  other  existing  lattice  Boltzmann  models  for  non-ideal  gases,  and  the 
similarities  and  differences  among  the  existing  models  are  explicitly  shown.  Sec.  6  concludes  the  paper.  A 
more  detailed  discussion  on  the  Enskog  equation  for  dense  gases  and  derivation  of  the  collision  term  leading 
to  non-ideal  gas  equation  of  state  are  provided  in  Appendix  A.  The  Chapman-Enskog  analysis  for  the  lattice 
Boltzmann  model  for  non-ideal  gases  is  demonstrated  in  Appendix  B.  In  Appendix  C,  the  forcing  term  in  the 
Boltzmann  equation  is  directly  derived  from  an  equilibrium  distribution  function  shifted  by  the  acceleration 
due  to  an  external  field.  This  derivation  provides  a  simple  and  clear  derivation  of  the  models  utilizing  the 
external  force  to  mimic  the  non-ideal  gas  effect. 

2.  Enskog  Equation  for  Dense  Gases.  The  Enskog  equation  [5,  30,  23]  explicitly  includes  the  radius 
of  colliding  particles,  ro,  in  the  collision  integral: 

(2.1a)  &/  +  £-V/  +  a-Vc/  =  J, 

(2.1b)  J  =  Jdiii  [g{x  +  r0r)f(x,  £')f(x  +  2 r0r,  £()  -  g(x  -  r0r)f(x,  £)f(x  -  2 r0f,  &)]  , 

where  /  is  the  single  particle  (mass)  distribution  function,  £  and  a  are  particle  velocity  and  acceleration, 
g  is  the  radial  distribution  function,  r  is  the  unit  vector  in  the  direction  from  the  center  of  the  second 
particle  of  f(x,  £1)  to  the  center  of  the  first  particle  of  f(x,  £)  at  the  instant  of  contact  during  a  collision, 
and  fi\  is  the  collisional  space  of  the  second  particle  of  f(x,  £1).  If  we  expand  the  collision  operator  J  in 
a  Taylor  series  about  x ,  use  the  BGK  approximation  [2,  23,  63,  40],  and  assume  the  fluid  to  be  isothermal 
and  incompressible,  we  obtain  the  following  equation  (details  refer  to  Appendix  A): 

(2.2a)  dtf  +  £•  V/  +  a- Vc/  =  -|[/  -  /<0)]  +  J' , 

(2.2b)  J'  =  -fw  bp  g(£,~  u)  -V  In  (p2g) , 

where  A  is  the  relaxation  time  which  characterizes  a  typical  collision  process,  and  /(0)  is  the  Maxwell  local 
equilibrium  distribution  function  [47]  given  by 

(2.3)  f(0)(p,  u,  6)  =  ^p0{2ne)~D/2  exp  [-(£  -  uf  /26  -  U(x)/9\  , 

where  D  is  the  dimension  of  the  momentum  space  p,  u ,  and  0  =  are  mass  density,  macro¬ 

scopic  velocity,  and  the  normalized  temperature,  respectively;  &#,  T,  and  m  are  the  Boltzmann  constant, 
temperature,  and  molecular  mass,  respectively;  U{x)  is  a  mean-field  external  potential  (per  unit  mass), 

(2.4)  po  =  ±JdxdSf 0) 
is  the  average  mass  density  in  the  system  of  volume  V,  and 

(2.5)  z(0)  =  f  dx  exp [—U(x)/0\ . 

V  Jv 

The  additional  collision  term,  J'  in  Eqs.  (2.2),  includes  the  volume  exclusion  effect  (see  Appendix  A  for 
details).  In  the  original  work  of  Enskog,  g  =  g(bp),  and  b  is  the  second  virial  coefficient  in  the  virial 
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expansion  of  the  equation  of  state  for  the  hard-sphere  system  [5].  The  hydrodynamic  moments,  i.e .,  mass 
density  p,  velocity  u ,  normalized  temperature  0 ,  and  energy  density  e  can  be  defined  as  follows: 

(2.6a)  p  =  Jd£f(0)  =  jd£f  =  ^p0exp(-U/Q) , 

(2.6b)  pu  =  jd&r  =  jdttf , 

(2-6c)  ^ pO  =  Jd£^(£-  u)2f(0)  =  Jd£^(£-  u)2f , 

(2.6d)  pe  =  Jd£  ^(€~u)2  +  U(x)  /(0)  =  jd£  ^(£-u)2  +  U(x)  f. 

The  acceleration  a  is  purely  due  to  an  external  field,  U(x),  i.e., 

(2.7)  a  =  £  =  -VU. 

For  the  Enskog  equation  or  the  revised  Enskog  equation,  both  global  [56]  and  local  [44,  45,  52]  iL-theorem 
can  be  proved. 

To  emphasize  that  the  acceleration  a  is  strictly  due  to  an  external  field  or  its  equivalent,  the  derivation  of 
the  Boltzmann  equation  via  the  BBGKY  hierarchy  is  briefly  reviewed  in  the  following.  Given  a  Hamiltonian 
system  of  N  particles,  the  7V-particle  distribution  function  Jn(x^n\  t)  satisfies  the  following  Liouville 

equation: 

N 

(2.8)  dtfN  +  fe* Vj  +  In  —  o, 

where  Vi  =  VXi,  =  {x\,  #2,  •  •  • ,  &n}  (similar  for  £( N )),  and 

(2.9)  &  =  a<  -  Y  VjVdlxi  -  xj\\)  =  -ViU  > 

and  di  =  —ViU  is  the  acceleration  due  to  external  force  (one-body  interaction  U),  and  Vij  =  K(||xi  —  Xj\\) 
is  the  two-body  interaction  potential  among  the  particles  (Vu  =  0).  By  distinguishing  the  external  field  and 
interparticle  interactions,  the  Liouville  equation  can  be  rewritten  as  the  following: 

N  N 

(2.10)  dtfN  +  ^  (£i* Vi  +  CLi-V^i)  f N  =  ^  • 

i=l  i= 1  jz£i 

The  reduced  distribution  function,  /m  for  1  <  M  <  N, 

fM  =  Jdx<N~MWN~M)fN, 

where  x^N~M^  =  {xn-m+u  &N-M+2,  •  •  • ,  xn}  (and  similar  for  £( N~M ))?  satisfies  the  following  equation 
[22,  23]: 

M  M  r 

(2.11)  dtfM  +  ^  +  ar%)  /m  =  (N  -  M)  ^  /  dXM+ld^M+lV M+\Vi,M+l  '  %/m+1  , 

i= 1  i= 1  ' 

The  above  equation  is  the  celebrated  BBGKY  hierarchy.  The  first  equation  in  the  hierarchy  for  the  single 
particle  distribution  function  f(x,  £,  t)  =  fi{x\,  £1,  t)  is: 

(2.12)  dtf  +  £•  V/  +  o-VJ  =  (AT  -  1)  J dx2d$2  V2U2  •  V€/2(a;,  £,  x2,  £2,  t) , 
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where  a  =  ai  =  —  Vi?7  is  purely  an  effect  due  to  the  external  field  U .  The  two-particle  distribution  function 
/2  can  be  decomposed  symbolically  as  the  following: 

(2.13)  f2(x i,  £1,  x2,  6,  t)  =  C2(xi,  x2,  t)f(x i,  £i,  t)f(x2,  £2,  t) , 

where  C^ix i,  #2,  t)  is  a  pair  correlation  function  [23,  35].  Depending  on  the  approximation  applied  to  the 
pair  correlation  C^ix  1,  #2,  t),  different  closure  would  result  from  the  kinetic  equation  (2.12)  [23,  35].  For 
example,  with  the  approximation  that  either  g2  =  1,  i.e ., 

f2(xi,  £1,  x2,  6,  t)  =  f(x  1,  £1,  i)/(a?2,  6,  0 , 

or  C^ixii  a?2,  £)  =  g(||#i  —  ®2||)j  where  g  is  the  radius  distribution  function  for  the  hard  sphere  system,  i.e., 
f2(x  1,  £1,  X2,  £2,  i)  =  3(11*1  -  *2||)/(*1,  £1,  t)f(x 2,  £2,  4 , 

the  first  equation  of  the  BBGKY  hierarchy,  Eq.  (2.12),  leads  to  either  the  Boltzmann  equation  or  the  Enskog 
equation,  respectively  [23,  35]. 

It  is  apparent  that  the  acceleration  a  is  due  to  a  self-consistent  external  field  which  is  one-body  interaction 
in  nature,  as  clearly  and  explicitly  illustrated  in  the  derivation  of  the  Boltzmann  equation  from  Liouville 
equation  via  BBGKY  hierarchy  [23,  30,  35,  36,  40].  In  other  words,  the  potential  U  in  the  Maxwellian 
defined  by  Eq.  (2.3)  only  represents  an  external  field  of  body-force  type,  and  this  self-consistent  mean-field 
interaction  should  not  be  confused  with  genuine  multi-body  interactions  which  take  place  in  non-ideal  gases. 
In  the  Boltzmann  equation,  all  the  interactions  among  particles  (multi-body  interactions)  involved  in  a 
collision  process  are  considered  in  the  collision  operator,  represented  by  collision  cross  section.  In  particular, 
the  collision  operator  reduces  to  a  parameter  A  of  the  single  relaxation  time  in  the  case  of  the  BGK  equation. 
In  the  limit  of  BGL,  the  interactions  among  particles  have  no  effect  other  than  changing  the  numerical  value 
of  the  viscosity,  as  clearly  illustrated  by  the  BGK  model.  Therefore,  non-ideal  gas  effects  are  not  included  in 
the  Boltzmann  equation.  To  exhibit  non-ideal  gas  effects  in  the  thermodynamical  limit,  the  finite  range  of 
the  interactions  among  particles  in  the  same  limit  (the  finite  size  effect  or  the  volume  exclusion  effect),  which 
causes  non-ideal  gas  effects,  must  be  explicitly  considered.  As  is  shown  in  detail  later  in  Sec.  5,  the  existing 
LBE  models  use  some  form  of  one-body  interaction  to  mimic  non-ideal  gas  effects.  While  this  approximation 
of  multi-particle  interaction  by  a  mean-field  one-body  interaction  seems  to  allow  LBE  models  to  simulate 
isothermal  non-ideal  gases  because  the  effect  of  pressure  and  forcing  can  be  distinguishable  in  the  momentum 
equation.  However,  this  is  no  longer  true  in  the  energy  equation.  Pressure  and  forcing  act  quite  differently 
in  the  energy  equation:  the  former  affects  the  energy  transport  as  PV  •  u  whereas  the  latter  does  work  as 
pa  u.  In  addition,  the  multi-particle  interactions  affect  the  heat  conductivity  whereas  the  forcing  does  not. 
This  suggests  some  inevitable  adverse  consequences  of  the  approximation.  The  only  way  to  correctly  model 
non-ideal  gases  is  to  at  least  include  the  “finite  size  effect”  explicitly.  And  the  Enskog  equation  is  one  such 
model. 

A  formal  solution  of  the  Enskog  equation  with  BGK  approximation,  Eq.  (2.2),  can  be  obtained  by 
integrating  along  characteristic  line  £  over  a  time  interval  of  length  8t  [36] : 

f{x  +  £5t  +  ^a<5t2,  €  +  aSt,  t  +  st)  =  e~St9/x  f(x,  £,  t)  +  ■Le~$tg/\  j  et'g/\  £  +  at' ,  t  +  t')  dt' 

(2.14)  -\-e~dtg/^  f  et  g/x  jt (x(P),  £  +  at' ,  t  +  t')  dt'  —  e~5t9^xa  •  f  el  9^x  V^f(x(tf),  £  +  at' ,  t  +  t')  dt' , 

Jo  Jo 

where  x(t')  =  x(t)  +  £t'  +  \<xt'2  is  the  (approximated)  trajectory  under  the  influence  of  the  external  field.1 
!We  neglected  the  term  of  t'2  in  [43]  because  it  does  not  affect  the  final  result. 

5 


The  approximation  is  made  by  the  assumption  that  the  acceleration  a  is  a  constant  locally.  Note  that  the 
above  equation  is  implicit  because  of  not  only  the  term  \^/,  but  also  the  time-dependence  of  hydrodynamic 
moments  p,  u ,  and  0  in  the  equilibrium  /(0),  and  J' . 

Our  derivation  of  the  lattice  Boltzmann  equation  is  based  upon  the  discretization  of  the  above  integral 
solution  of  the  Enskog  equation.  In  what  follows,  we  show  that  the  lattice  Boltzmann  equation  is  an  explicit 
finite  difference  scheme  for  solving  the  above  integral  solution  of  the  Enskog  equation. 

3.  Derivation  of  Lattice  Boltzmann  Equation. 

3.1.  Discretization  in  time.  By  using  the  mean- value  theorem,  we  can  rewrite  the  integral  solution 
of  the  BGK  Enskog  equation,  Eq.  (2.14),  as  follows: 

/(*  +  £6t  +  \a82,  £  +  aSt,  t  +  St)  =  e~^/x  f(x,  £,  t)  +  je-s^°/xf  °\x(e8t),  £  +  aeSt,  t  +  eSt)  8t 
(3.1)  +e~St{1~e)9/x  J' {x(e5t),  £  +  aeSt,  t  +  eSt)  St  -  e~St(-1~e')9/xa-Vif(x(eSt),  £  +  aeSt,  t  +  eSt)  8t , 


where  e  is  a  constant  between  0  and  1,  and  x(eSt)  =  x(t)  +  $,eSt  +  \a{e8t)2 .  If  we  assume  that  St  is  small 
enough  and  /(0),  J',  and  are  smooth  enough  locally  in  phase  space,  we  can  neglect  the  terms  of  the 
order  0(82)  or  smaller  in  the  Taylor  expansion  of  Eq.  (3.1),  and  obtain 

(3.2)  f(x  +  £8t,  4,  t  +  St)~  f(x,  4,  t)  =  -~[f(x,  £,  t)  -  f  0)(x,  £,  t)]  +  J'(x,  t)  5t  -  a-^f(x,  t )  5t , 

T 

where  r  =  X/8t  is  the  dimensionless  relaxation  time.  It  is  obvious  that  the  accuracy  of  the  above  equation 
is  only  first  order  in  time.  Consequently  the  accuracy  of  the  lattice  Boltzmann  equations  derived  from  the 
above  equation  is  also  first  order  in  time  in  principle. 


3.2.  Low  Mach  number  expansion  and  phase  space  discretization.  There  are  two  steps  in  the 
derivation  of  lattice  Boltzmann  equation  from  Eq.  (3.2):  (a)  construction  of  an  appropriate  equilibrium  dis¬ 
tribution  function,  and  (b)  a  coherent  discretization  of  phase  space.  For  the  isothermal  case,  the  equilibrium 
distribution  function  can  be  obtained  by  truncation  of  the  Taylor  expansion  of  /(0)  up  to  the  second  order 
in  u : 


(3.3) 


/(eq)  =  ^-p0(2n9)-D/2  exp(—U/0)  exp(-£2/2<?) 
= 


\  |  (*•«)  |  (g-«)2 
6  202  20 


(^w)  (£-u)2 

0  202 


u 

20 


where 


(3.4)  o>(£)  =  (2n0)  13/2  exp  (-£2/2i 9). 

The  phase  space  discretization  has  to  be  done  in  such  way  that  not  only  all  the  hydrodynamic  moments, 
but  also  their  corresponding  fluxes  are  preserved  exactly.  This  implies  that  the  following  quadrature  must 
be  evaluated  exactly: 

(3.5)  Jd£€kf{eq),  0  <  fc  <  3 , 

for  isothermal  models.  (Here  we  require  that  not  only  all  the  hydrodynamic  moments,  but  also  the  corre¬ 
sponding  fluxes,  are  computed  exactly  by  the  quadrature.  This  requirement  is  perhaps  more  stringent  than 
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necessary  because  energy  flux  is  usually  not  considered  in  the  isothermal  case.)  Because  of  the  second  order 
polynomial  contained  in  / (eq)  given  by  Eq.  (3.3),  the  quadrature  which  must  be  evaluated  exactly  is: 

(3.6)  Jd££kex p(-£2/20) ,  0  <  k  <  5  . 

Because  of  the  exponential  function  in  the  above  integral,  Gaussian  quadrature  [7]  is  a  natural  choice  for 
the  evaluation  of  the  integral.  With  a  k- th  order  polynomial  ^k{x)  of  x,  Gaussian  quadrature  defined  by  the 
following  equation, 

/oo  n 

dx  ipk  (x)  e_x2/2  =  ^2  Waxpk  (xa) , 

-oo  1 


is  exact  for  0  <  k  <  2n  —  1,  where  Wa  and  xa  are  the  weights  and  the  abscissas  of  the  quadrature,  respectively. 

3.3.  The  forcing  term.  The  forcing  term,  a  -\^/,  must  be  constructed  explicitly  in  the  lattice  Boltz¬ 
mann  equation.  We  use  the  moment  constraint  to  construct  this  term.  The  moments  (up  to  the  second 
order)  of  the  forcing  term  are: 

(3.8a)  Jd^a-^f  =  Jd$a-\7^f0)  =  0, 

(3.8b)  Jdtta-\7€f  =  Jdtta-W»  =  -pa, 

(3.8c)  Jd£ a-^f  =  Jd4 a- V€/(0)  =  -p(aiUj  +  ajUi) . 

Here,  we  have  noted  that  /  can  be  replaced  (or  approximated)  by  /(0)  without  affecting  the  moments  of  the 
forcing  term  up  to  the  second  order  in  £  —  in  general  the  replacement  of  /  by  /(0)  does  not  hold  for  the 
moments  higher  than  the  third  order  in  £.  This  is  owing  to  the  fact  that  /  and  /(0)  have  exactly  the  same 
conserved  (or  hydrodynamic)  moments,  a  constraint  on  the  normal  solution  of  the  Boltzmann  equation  in 
the  Chapman-Enskog  analysis. 

The  forcing  term,  can  be  written  in  terms  of  an  expansion  in  £  as  follows: 

(3.9)  a- W  =  pu(£)  [c<°>  +  cf ]ti  +  +  ■  •  •]  , 


where  the  Einstein  notation  of  summation  for  the  repeated  Roman  indices  i,  j,  . . . ,  is  used.  The  first  few 
coefficient  -  can  eas^Y  obtained  by  using  the  moment  constraints  given  by  Eqs.  (3.8)  if  the  above 
expansion  is  truncated.  With  the  truncated  expansion  up  to  the  second  order  in  £  and  the  first  order  in  u , 
we  obtain 


(3.10a) 

(3.10b) 

(3.10c) 


,(o)  = 


f2 

ST 


(1)  1 
Ci  =-72a*’ 
vr 

cif  =  -X7T  (aiu3  +  aiui)  > 

ZS  T 


where  £ x  =  \/0  is  proportional  to  the  thermal  velocity  of  a  particle  at  temperature  T.  Therefore,  up  to  the 
order  of  0(u)  and  0(£2),  we  have 


(3.11) 


«-V€/  =  -pto(€)£T2  [(£-u)+tT2(£-u)t,]-a. 
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Note  that  in  the  above  expansion,  only  the  terms  up  to  the  first  order  in  u  have  been  retained,  because  there 
is  an  overall  factor  of  St  in  the  forcing  term,  as  indicated  in  Eq.  (3.2).  And  St  is  of  0{u)  in  the  Chapman- 
Enskog  analysis  for  the  lattice  Boltzmann  equation  (see  Appendix  B  for  explanation).  It  should  be  stressed 
that  every  term  in  the  Boltzmann  equation  must  be  treated  equally  in  terms  of  maintaining  the  accuracy. 
Specifically  speaking,  the  expansion  of  the  forcing  term  must  be  of  the  second  order  in  £  and  the  same  in  the 
small  expansion  parameter  Stj  in  order  to  be  consistent  with  the  expansion  of  the  equilibrium.  It  should  be 
noted  that  there  are  other  methods  to  compute  the  expansion  of  the  forcing  term.  Up  to  the  second  order  in 
and  the  first  order  in  u ,  the  expansion  of  a-\^/(0)  is  identical  to  that  of  a-\^/  because  of  the  constraints 
given  by  Eqs.  (3.8).  Therefore,  the  result  of  Eq.  (3.11)  can  be  obtained  by  computing  a-\^/(0)  explicitly. 

It  should  be  pointed  out  that  there  is  another  way  to  include  the  effect  of  the  forcing  due  to  an  external 
field.  Assuming  that  multi-body  interactions  among  the  particles  in  the  system  are  of  short-range  and  the 
mean-free-path  of  a  particle  is  much  larger  than  the  interaction  range,  then  a  particle  is  accelerated  only  by 
the  external  field  between  collisions.  Thus  the  net  effect  of  the  acceleration  due  to  the  external  field  during 
the  mean-free-time  is  an  increment  of  particle  velocity.  Therefore  one  can  use  an  equilibrium  distribution 
function  with  a  velocity  shift  to  account  for  the  effect  of  the  forcing  due  to  the  external  field  [66],  i.e., 
/(0)(p,  u ,  6)  becomes  /(0)(p,  u  —  arSt ,  6)  at  the  presence  of  the  external  field.  Naturally,  the  accelerated 
equilibrium  distribution  function  /(0)(p,  u  —  arSt ,  0)  leads  to  a  forcing  term  in  the  lattice  Boltzmann  equation 
when  discretized  (see  Appendix  C  for  details).  It  should  be  noted  that  these  two  approaches  are  equivalent 
up  to  the  first  order  in  8t.  At  higher  order  of  St  the  velocity  shift  in  the  equilibrium  distribution  will  introduce 
nonlinear  terms  which  are  different  from  what  derived  from  the  continuous  equation. 

3.4.  2D  nine-velocity  model  on  square  lattice.  We  now  use  the  2D  nine-velocity  LBE  model  on 
a  square  lattice  space  as  a  concrete  example  to  illustrate  our  discretization  scheme.  A  Cartesian  coordinate 
system  (in  £-space)  is  used  in  this  case,  and  accordingly  we  set  t/>(£)  =  Thus  the  quadrature  needed 

to  be  evaluated  is  the  following: 


(3.12) 


-  ck+l+2hh 


/  =  £ 


where 

(3.13) 


h 


e-<2/2, 


and  (  =  £x  / or  £v/£t-  Naturally,  the  third-order  Hermite  formula  [7]  is  the  optimal  choice  to  evaluate 
Ik  for  the  purpose  of  deriving  the  9-bit  LBE  model,  i.e.,  Ik  =  Y^j= i  u j(j  •  The  three  abscissas  (Q)  and  the 
corresponding  weights  (ouj)  of  the  quadrature  are: 


(3.14a)  Ci  = -a/3,  C2=0,  C3  =  a/3, 

(3.14b)  uji  =  VT/6  ,  lJ‘2  =  2^7r/3  ,  cj3  =  y/ir/Q  . 

Then,  the  integral  of  Eq.  (3.7)  becomes: 

4  8 

^2^(0)  +  ^  +  E  ‘‘•’lVKCa)  , 

ol— 1  a=5 

where  is  the  zero  velocity  vector  for  a  =  0,  one  of  the  vectors  of  y/S  (±1,  0)  and  y/S  (0,  ±1)  for 
a  =  1-4,  and  one  of  the  vectors  of  y/S  (±1,  ±1)  for  a  =  5-8.  Note  that  the  above  quadrature  is  exact  for 
the  integral  defined  by  Eq.  (3.13)  when  k  <  5. 


(3.15) 


1  =  2$ 
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Now  momentum  space  is  discretized  with  nine  discrete  velocities  =  0,  1,  •  •  • ,  8}.  To  obtain  the 

9-bit  model,  configuration  space  is  discretized  accordingly,  i.e.,  it  is  discretized  into  a  square  lattice  space 
with  a  lattice  constant  Sx  =  y/3  St .  It  should  be  stressed  that  the  temperature  T  (or  0)  is  a  constant  here 
because  we  are  only  dealing  with  an  isothermal  model.  We  can  therefore  choose  Sx  to  be  a  fundamental 
quantity  instead,  thus  y/S  =  c  =  Sx/St,  or  0  =  =  c2/ 3.  Thus,  phase  space  is  discretized  coherently: 

the  discretizations  of  the  velocity  space  and  the  configuration  space  are  closely  coupled  together.  This  is  one 
feature  of  the  lattice  Boltzmann  equation  distinctive  from  other  finite  difference  schemes. 

By  comparing  Eqs.  (3.7)  and  (3.15),  we  can  identify  the  weights  defined  in  Eq.  (3.7): 

(3.16)  Wa  =  2ir  &  exp(£a/2 &)  > 


where 


(3.17) 


wn  =  < 


4/9, 

1/9, 


a  =  0, 

a  =  1,2,  3,  4, 


[  1/36,  a  =  5,  6,  7,  8. 
Then,  the  equilibrium  distribution  function  for  the  9-bit  model  is: 

fieq)  =Waf^(x,£a,t)=wap 


(3.18) 
where 

(3.19) 


1  3(ea  ■  u)  9(ea  ■  u )2 


2c4 


3u2 

2c2” 


en  =  < 


f  (0,0), 

(cos  (j)a,  sin  (j)a)  c, 


a  =  0 , 
a  =  1  2,  3,  4, 
a  =  5,  6,  7,  8, 


[  (cos  </>a,  sin  (f>a)  V2c, 

and  (J)a  =  (a  —  l)7r/2  for  a  =  1-4,  and  (J)a  =  (a  —  5)7r/2  +  7r/4  for  a  =  5-8. 

3.5.  Discretized  forcing  term.  Applying  the  same  discretization  to  the  forcing  term  of  Eq.  (3.11), 
we  have  the  discretized  forcing  term  for  the  nine- velocity  model: 


(3.20) 


Fa  =  -3  wap 


~j(ea  —  U)  +  3 


(ea-u) 


[c  c* 

The  forcing  in  the  above  equation  satisfies  the  following  constraints: 


•a . 


(3.21a) 

O 

II 

W' 

(3.21b) 

^  ^  Fa  —  p  CL  , 

OL 

(3.21c) 

^  ^  Ca,i&a,j  Fa  = 

The  above  constraints  are  the  discrete  counterpart  of  Eqs.  (3.8).  If  only  the  first  two  moment  equations  in 
Eqs.  (3.21)  are  to  be  satisfied,  and  the  third  constraint  of  Eq.  (3.21c)  is  replaced  by 


(3.22) 

the  forcing  term  thus  reduces  to 

(3.23) 


^  ^  Ca,i^a,j  Fa  —  0  , 


■jji  _  9„..  „{ea-a) 

Fa  —  5  Wa  p  2 
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The  above  forcing  term  is  what  has  been  often  used  for  constant  body  force  in  the  literature  [41,  42].  One 
adverse  consequence  of  using  the  above  forcing  term  is  that  the  Galilean  invariance  is  lost  if  a  is  not  a 
constant  in  space.  In  addition,  the  work  by  forcing,  pa  •  u,  does  not  appear  in  the  energy  balance  equation, 
and  thus  leads  to  an  incorrect  energy  balance  equation.  As  shown  in  Sec.  5,  the  forcing  terms  of  similar 
forms  are  used  to  produce  various  non-ideal  gas  effects  in  previous  models  [59,  60,  62,  61,  51,  50,  27]. 

3.6.  The  lattice  Boltzmann  equation.  The  additional  collision  term,  J '  of  Eq.  (2.2b),  can  be  ex¬ 
plicitly  written  in  the  discrete  form: 

(3.24)  J'a  =  — /<0)  bpg(ea  —  u)-'V\n(p2g) . 

With  the  discretized  J'  included,  the  lattice  Boltzmann  equation  for  dense  gases  is: 

fa(x  +  ea8t,  t  +  8t)  —  fa(x,  t)  =  [/<*(»,  t)  -  fieq)(x,  £)] 

(3.25)  -bpgf^(x,  t)  (ea  -  u)-V(p2g)  -Fa8t, 

where  the  forcing  term,  Fa,  is  given  by  Eq.  (3.20).  The  hydrodynamic  moments  in  the  lattice  Boltzmann 
models  are  given  by 

(3.26a)  p  =  ^/a  =  ^/(eq), 

a  a 

(3.26b)  pu  =  ^  eafa  =  y^ea/^eq) , 

a  a 

(3.26c)  p6=\  _  =  \  _  u)2fieq)  ■ 

a  a 

The  additional  collision  term,  J^,  involves  the  density  gradient,  Vp,  which  can  be  explicitly  computed  by 
either  the  second  order  central  differencing 

ea-Vp(x)8t  =  ^\p{x  +  eaSt)  -  p(x  -  eaSt)] , 

or  the  first  order  differencing 


ea-Vp{x)  St  =  p(x  +  eaSt)  -  p(x) . 


The  other  alternative  would  be  to  construct  a  collision  term  similar  to  the  original  Enskog  collision  term 
given  by  Eq.  (2.1b)  without  the  Taylor  expansion  in  space. 

4.  The  Hydrodynamics  and  Thermodynamics.  Through  the  Chapman-Enskog  analysis  (see  Ap¬ 
pendix  B  for  the  details),  the  hydrodynamic  equations  of  the  lattice  Boltzmann  model  for  dense  gases,  given 
by  Eq.  (3.25)  with  the  equilibrium  of  Eq.  (3.18),  are: 


(4.1a) 

(4.1b) 

where  the  viscosity 


(4.2) 


dtp  +  V-(pu)  =  0, 

dtu  +  u-Vu  =  — -VP  +  vV2u  +  a  , 
P 


v  =  — - - cdx  . 

6  g 
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and  the  pressure  (or  the  equation  of  state)  is  given  by 
(4.3)  P  =  p6(l  +  bpg). 


With  the  above  equation  of  state,  the  sound  speed,  cs,  becomes 


(4.4)  cl=0 

For  ideal  gas,  b  =  0  and  g  =  1,  P,  z/,  and  cs  recover  the  previous  results  for  ideal  gas.  The  dependence  of 
the  viscosity,  z/,  on  g  can  be  removed  by  replacing  g  in  the  BGK  collision  term  by  1. 

Although  in  the  original  work  of  Enskog  [5],  g  only  accommodates  the  volume  exclusion  effect,  or  repulsive 
interaction,  in  the  gas  of  hard  spheres,  however  there  is  no  reason  to  prohibit  the  inclusion  of  a  more  general 
interaction.  Indeed,  g  can  be  somewhat  arbitrary,  depending  on  the  interaction.  The  radial  distribution 
function  g  provides  a  freedom  to  alter  the  transport  coefficients  (z/  and  cs)  as  well  as  the  equation  of  state. 
However,  it  should  be  stressed  that  there  are  bounds  to  this  freedom.  From  Eq.  (4.2),  it  becomes  obvious 
that  the  model  is  stable  if  and  only  if  r  >  g/ 2.  This  suggests  that  g  also  affects  the  numerical  stability 
of  the  system.  In  addition,  the  sound  speed  can  be  changed  by  g.  But  one  must  not  expect  to  achieve 
cs  >  c  =  Sx/St  or  the  basic  principle  of  physics  would  be  violated,  because  c  limits  the  speed  of  information 
propagation  in  the  LBE  system.  Therefore,  there  are  bounds  to  the  values  of  g  and  derivative  of  p2g. 

With  the  equation  of  the  state  given,  the  Helmholtz  free  energy  density  can  be  given  by: 


1  +  *E<*> 


(4.5) 


Mp) 


=  pJ  ~^dp  =  p0  In p  +  b  J  gdp 


because 


(4.6)  F  = 

And  the  radial  distribution  function  g  can  also  be  computed  from  either  P  or  ip.  That  is,  with  either  P  or  ip 
given,  one  can  derive  all  the  relevant  thermodynamic  quantities  from  the  free  energy  ip.  For  example,  given 
the  van  der  Waals  equation  of  state: 


(4.7) 


P  =  p  0 


1 


[(1-bp)  0P 


where  parameter  a  accounts  for  the  mean  result  of  attractive  potential  among  particles  [5] ,  and  according  to 
Eqs.  (4.3)  and  (4.7)  the  radial  distribution  function  g  is 

1  a 

(1  —  bp)  bO 


9  = 


And  the  corresponding  free  energy  density  is 


ip  =  pO 


In 


~0P 


A~bP, 

With  the  free  energy  and  the  equation  of  state  defined,  the  Maxwell  construction  [32]  to  determine  the 
co-existence  curve  becomes  physically  meaningful  and  consistent.  Nevertheless,  care  must  be  taken  in  con¬ 
ducting  the  Maxwell  construction  in  the  discretized  situation.  The  phenomena  of  liquid-gas  phase  transition 
can  be  simulated  by  the  model  by  charging  the  value  of  5^ 9&P  (or  simply  just  b)  in  the  free  energy  density 


ip  relative  to  the  temperature  #,  as  indicated  by  Eq.  (4.5).  Bear  in  mind  that  the  temperature  0  cannot 
be  changed,  because  it  is  a  fixed  constant  in  the  isothermal  LBE  models.  It  should  be  noted  that  g  should 
be  computed  with  a  given  potential  in  principle.  The  above  manipulation  to  obtain  g  is  not  based  upon 
principles  of  physics.  Also,  the  use  of  the  free  energy  adds  nothing  to  the  physics  of  the  model,  it  only  reflects 
a  matter  of  custom  or  preference. 
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5.  Analysis  of  Some  Existing  Models.  What  we  propose  in  this  work  is  a  systematic  construction 
of  the  lattice  Boltzmann  equation  in  a  consistent  and  a  priori  fashion,  with  the  premise  that  the  continuous 
Boltzmann  equation  is  adequate  to  describe  underlying  physics  of  the  systems  of  interest.  In  particular,  for 
non-deal  gases,  one  must  use  the  Enskog  equation  for  dense  gases  instead  of  the  original  Boltzmann  equation 
for  dilute  gases.  In  light  of  this  view  point,  a  survey  of  the  existing  LBE  models  for  non-ideal  gases  is  now 
in  order.  We  discuss  two  lattice  Boltzmann  models  for  non-ideal  gases  which  are  independently  proposed  by 
Shan  and  Chen  [59,  60],  and  by  Swift,  Osborn,  and  Yeomans  [62,  61].  In  spite  of  the  significant  differences  in 
their  appearances  and  in  technical  details,  these  models  do  share  one  common  feature  in  their  constructions 
of  the  lattice  Boltzmann  model  for  non-ideal  gases:  The  derivation  of  the  lattice  Boltzmann  models  is  mainly 
accomplished  by  constructing  a  phenomenological  equilibrium  distribution  function  which  can  accommodate 
non-ideal  gas  effects  and  satisfies  all  the  conservation  constraints,  and  therefore  leads  to  hydrodynamics.  In 
what  follows,  we  shall  analyze  these  two  models  and  explicitly  demonstrate  the  difference  between  the  model 
derived  in  this  paper  and  the  aforementioned  two. 

5.1.  Model  with  interacting  potential.  In  the  model  proposed  by  Shan  and  Chen  [59,  60],  a  local 
density-dependent  potential  U(p(x))  ~  Q6ip2(p)  is  explicitly  given,  where  Q  is  the  interaction  strength  and 
pj  is  an  arbitrary  function  of  density  p.  The  change  of  the  particle  velocity  £  (not  the  macroscopic  velocity 
u)  due  to  U(x)  is 


5£  =  —VU(x)  rSt  =  a  rSt , 


and  Su  =  —5£  is  explicitly  substituted  into  the  equilibrium  distribution  function,  i.e., 


f^=wap  | 


1  + 


3  [ea-(u  —  ardt)\  9  [ea  ■  (u  —  arSt)]2  3  (u  —  adt)2 


2c4 


(5.1) 


=  wap 


1  +  3  (ea-u)  +  9(ea-u)2 


2c4 


3  u2 

2^ 


-3  wap 


2  c2 


—  (ea  —  u)  +  3 


(ea-u) 


■a  r8t 


W*P 


a_  _  ( ea  a )2 

r-2  rA 


_2  c 2 
T  0+  . 


In  the  above  result,  the  first  part  is  the  usual  equilibrium  distribution  function  which  has  an  ideal  gas  equation 
of  state  built  in,  and  the  second  part  accounts  for  the  interaction  or  non-ideal  gas  effects,  which  is  identical 
to  the  forcing  term  given  by  Eq.  (3.20),  produced  by  the  forcing  term  a-V^f  in  the  streaming  operator. 
By  combining  the  forcing  with  the  pressure  term  in  the  Navier-Stokes  equation,  the  equation  of  the  state 
becomes  P  =  p  [0  +  U{p)\.  Thus,  non-ideal  gas  effects  are  obtained  through  the  phenomenological  potential 
U.  To  achieve  the  purpose  of  mimicking  non-ideal  gas  effects,  the  leading  term  in  the  density  expansion  of  U 
has  to  be  of  second-order  in  p,  i.e .,  U  ~  Q6p2 ,  or  ip  ~  p,  as  specifically  indicated  in  [59,  60].  Obviously,  the 
potential  U{p{x))  is  intended  to  be  the  inter-particle  interaction.  However,  it  is  mathematically  implemented 
as  an  external  field  such  that  its  only  effect  is  producing  a  term  VU  in  the  momentum  equation  [59,  60].  The 
consequence  of  this  conceptual  confusion  is  that  the  energy  balance  equation  is  incorrect,  because  the  result 
of  an  external  field  is  the  work  of  pa  •  u,  while  the  result  of  the  inter-particle  interaction  is  a  heat-transfer 
due  to  the  viscous  effect,  as  shown  in  Appendix  B.  Specifically,  in  the  energy  equation,  the  correct  term 
related  to  the  pressure  is  PV  •  u ,  where  pressure  is  exactly  the  one  that  appears  in  the  momentum  equation. 
However,  with  the  one-body  interaction,  this  becomes  pO  V-  u  +  pVU  •  u ,  i.e.,  the  equation  of  state  is  not 
the  same  in  the  momentum  equation  and  the  energy  equation.  Furthermore,  the  third  part  in  Eq.  (5.1), 
which  is  proportional  to  S2  and  nonlinear  in  a,  is  omitted.  This  term  can  be  significant  when  St  is  set  to 
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unity,  as  it  is  a  common  practice  in  the  lattice  Boltzmann  simulations.  It  should  also  be  pointed  out  that 
the  viscosity  in  this  model  remains  intact  —  it  is  not  affected  by  the  potential  U . 

We  have  also  noted  a  recent  attempt  to  theoretically  justify  the  model  of  Shan  and  Chen  [27].  With 
some  crude  approximations  [27],  He  et  al.  showed  that  a  desirable  forcing  term  to  mimic  non-ideal  gas 
effects  is  Fa  ~  fa  (ea  —  u)-FSt,  where  F  ~  —  W  —  bpOg'Vln(p2g)J  and  V  =  —2 ap  —  kV2 p  accounts  for  the 
attractive  part  in  the  inter-particle  interaction.  Without  any  surprise,  this  model  reproduces  an  anticipated 
non-ideal  gas  equation  of  state,  P  «  pO  (1  +  bpg )  +  pH,  and  avoids  the  nonlinear  term  of  S2  in  the  model  of 
interacting  potential,  as  expected.  However,  the  energy  balance  equation  from  this  model  is  still  incorrect, 
due  to  the  similarity  to  the  previous  model.  It  should  also  be  noted  that  it  is  conceptually  incorrect  to  write 
the  pressure  as  P  «  pO  (1  +  bpg)  +  pV.  One  correct  way  to  generalize  the  van  der  Waals  equation  of  state 
is  writing  it  as  the  following  [4]: 

(5.2)  P  +  ap2  =9p(l  +  bpg), 


where  parameter  a  is  related  to  the  two-body  interaction  potential  by: 


(5.3) 


a  = 


3  m2  J  dr 


In  the  Enskog  equation,  g  is  obtained  from  a  hard-sphere  gas,  thus  the  attractive  potential  has  to  be  inserted 
through  parameter  a. 

It  is  clear  that  the  inter-particle  interactions  are  conceived  mathematically  as  external  fields  in  the  afore¬ 
mentioned  models.  Perhaps  the  only  plausible  justification  for  this  view  is  that  the  inter-particle  interaction 
can  be  approximated  by  a  self-consistent  one-body  interaction  field  (as  in  the  Vlasov  approximation  for 
Coulomb  gases).  In  this  case, 


(5.4)  J dx2d$, 2  V2V12  ■  V(f2(x,  £,  x2,  t)V(f(x,  t)  ■  J dx2d£2  f(x2,  t)V2V12  =  V(f  ■  VU , 

where  the  Boltzmann  approximation  has  been  invoked,  and 

v£7  =  [ dX2d^2  f(x 2,  £ 2 ,  f)V2Vl2 


defines  the  self-inconsistent  mean-field  potential  U.  This  approximation  is  justified  for  rarefied  collisionless 
plasma  with  Coulomb  interactions,  and  is  simply  inappropriate  for  non-ideal  gas  systems. 


5.2.  Model  with  free  energy.  A  comparison  with  the  model  proposed  in  [62,  61]  is  slightly  more 
elaborate.  Stressing  the  consistency  of  thermodynamics  in  the  lattice  Boltzmann  equation  and  being  inspired 
by  Cahn-Hilliard’s  model  for  surface  tension  [3,  13],  the  model  proposed  by  Swift  et  al.  [62,  61]  starts  with 
a  free  energy  functional: 

(5.5)  y  =  Jdx[^\\Vp\\2  +  ^(p)]  , 

where  is  the  bulk  free  energy  density.  The  free  energy  functional  in  turn  determines  the  diagonal  term  of 
the  pressure  tensor: 

(5-6)  p  =  p—-y=p-KpV2p--\\Vp\\2, 

where 

(5.7)  p  =  p^f--xl, 
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is  the  equation  of  state  of  the  fluid.  The  full  pressure  tensor  is  given  by 
(5.8)  P ij  =  P  Sij  +  k  dip  djp . 

With  the  pressure  tensor  given,  the  equilibrium  distribution  function  is  constructed  such  that  not  only  it 
satisfies  the  conservation  constraints  of  Eqs.  (3.26),  but  it  also  produces  the  above  pressure  tensor  by  enforcing 
additional  constrains  fa°^ea,iea,j  —  P  ij  [62,  61],  even  though  Eq.  (5.8)  is  fundamentally  incorrect. 

To  make  our  analysis  as  transparent  as  possible,  it  is  crucial  to  make  an  explicit  connection  between  the 
model  of  Swift  et  al.  to  the  model  of  Shan  and  Chen.  The  equilibrium  distribution  function  in  the  model  of 
Swift  et  al.  with  triangular  lattice  [62,  61]  can  be  rewritten  as  follows: 

1  +  (ea  ■  u)  +  2(ea  ■  u )2  -  ^u2 

(5-9)  +^{(e2a,x  -e2aJ  [(dxp)2  -  (dyp)2]  +2 ea,xea,ydxpdyp}  -  ^pV2p+^[pip'(p)  -  ip(p)  -  p] 

where  c  =  Sx/St  is  assumed  to  be  unity.  The  term  in  brackets  [  ]  is  nothing  but  the  usual  equilibrium 
distribution  function  of  the  7-bit  model  [25,  26].  The  term  in  bracket  {  }  is  an  expression  for  the  tensor 
Eij  =  ( eajdip)(eajdjp )  written  in  terms  of  a  traceless  and  an  off-diagonal  part  with  correct  symmetry 
such  that  all  the  terms  proportional  to  k  reduce  to  the  term  K[pV2p  +  ||Vp||2/2]  in  the  diagonal  part  of  the 
pressure  tensor,  given  by  Eq.  (5.6).  This  term  is  directly  taken  from  Cahn-Hilliard’s  model  and  it  induces 
surface  tension  due  to  density  gradient  in  addition  to  the  part  due  to  the  (non-ideal  gas)  equation  of  state, 
but  it  does  not  contribute  to  the  hydrodynamic  pressure  (or  the  equation  of  state).  The  non-ideal  gas  part 
in  the  equation  of  state  is  contained  in  the  last  part  the  above  equation:  [p'ip'(p)  —  'ip(p)  —  p\/ 3,  which  can 
be  written  in  a  density  expansion  in  general: 

(5.10)  $ip=  t[p-0'(p)  -V’(P)  ~  P]  =  0p2(B  +  Cp+  •••) 

where  coefficients  £?,  (7,  •  •  •  are  virial  coefficients.  We  have  noted  that  only  the  leading  term  in  the  density 
expansion  of  (p,  Bp2 ,  is  needed  in  order  to  capture  the  non-ideal  gas  effects,  for  this  term  not  only  leads  to 
a  non-ideal  gas  equation  of  state,  but  also  provides  all  the  necessary  terms  to  control  the  surface  tension 
in  Cahn-Hilliard’s  model,  because  VVp2  =  2 (pV2 p  +  ||Vp||2).  By  comparing  Eq.  (5.9)  to  Eq.  (5.1),  the 
connection  between  the  model  of  Swift  et  al.  and  the  model  of  Shan  and  Chen  becomes  obvious  if  Eq.  (3.23) 
for  the  forcing  term  is  used  and  the  equivalence  of  0  ip  =  ea  •  VU  ~  Fa  is  established.  Thus,  the  model  of 
Swift  et  al.  uses  Eq.  (3.23)  for  the  forcing  term,  which  is  only  valid  for  a  constant  body  force,  whereas  the 
model  of  Shan  and  Chen  uses  Eq.  (3.20)  for  the  forcing  term.  The  interaction  strength  in  the  model  of  Swift 
et  al.  is  proportional  to  temperature  0  (which  is  incorrect  for  isothermal  fluids),  whereas  in  the  model  of 
Shan  and  Chen,  it  is  proportional  to  a  constant  Q.  Since  the  connection  of  the  two  models  can  be  explicitly 
established,  all  the  analysis  in  the  previous  subsection  can  be  immediately  applied  to  the  present  model. 

It  should  be  pointed  out  that  at  the  level  of  the  Boltzmann  equation,  the  density  gradient  term,  ||  Vp||2,  in 
the  free  energy  functional  has  no  justification  whatsoever  within  the  framework  of  Chapman-Enskog  analysis. 
[In  fact,  the  density  gradient  Vp  can  only  appear  in  the  second  order  solution  of  /  (the  Burnett  equation 
[5])  in  the  Chapman-Enskog  analysis,  which  is  beyond  the  Navier-Stokes  equations.]  It  is  clear  that  the 
analysis  presented  in  [62,  61]  and  the  subsequent  work  [51,  50]  does  not  follow  the  Chapman-Enskog  analysis; 
therefore,  it  cannot  lead  to  a  mathematically  valid  derivation  of  macroscopic  equations  from  the  mesoscopic 
equation.  Consequently,  the  derivation  of  macroscopic  equations  from  the  model  becomes  dubious,  and  the 
model  suffers  from  a  number  of  unexpected  adverse  effects  which  are  discussed  next. 
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First,  the  model  lacks  Galilean  invariance,  mainly  because  of  the  forcing  term  Fa  ~  0  ip.  This  defect 
can  be  fixed  by  using  a  correct  forcing  term  according  to  Eq.  (3.20):  Fa  p[(ea  -u)+  3(ea  •  u)ea]  •  Vy>. 
However,  there  are  other  terms  involving  the  density  gradient  which  also  cause  the  problem  of  lack  of  Galilean 
invariance  [31].  Second,  the  ratio  between  the  number  of  rest  particles  and  the  number  of  moving  particles 
depends  on  the  local  density  gradient.  It  can  be  shown  that  this  ratio  is  related  to  temperature,  because  in 
the  two-speed  system,  the  width  of  the  equilibrium  distribution,  which  is  the  temperature,  is  determined  by 
this  ratio.  [To  be  exact,  according  to  the  definition  of  a;,  Eq.  (3.4),  the  ratio  u(0)/u(c)  =  ec2/261,  which  must 
be  a  constant  in  isothermal  fluids.]  This  means  that  the  temperature  may  vary  locally  depending  on  the 
density  gradient  while  the  model  is  claimed  to  simulate  an  isothermal  fluid.  Again,  this  problem  can  be  fixed 
by  using  the  correct  forcing  term  aforementioned.  In  addition,  the  model  cannot  lead  to  the  correct  energy 
balance  equation  for  the  very  same  reason,  as  shown  in  the  previous  subsection,  that  the  terms  related  to  the 
free  energy  can  be  considered  as  a  body-force  due  to  a  thermodynamic  potential  (free  energy  in  the  model) 
which  is  a  mean-field  quantity. 

We  stress  that  the  conceptual  difficulty  caused  by  the  confusion  between  body-force  and  interaction 
cannot  be  circumvented  by  technical  tricks  such  as  using  a  correct  forcing  term,  or  including  higher  order 
terms  in  u  in  the  equilibrium  distribution  function  faq\  The  reason  is  obvious:  as  long  as  a  mean-field 
potential,  whether  an  interaction  or  a  free  energy,  is  employed  to  mimic  non-ideal  gas  effects,  the  constraints 
of  Eqs.  (3.8)  must  be  satisfied,  therefore  the  inconsistency  of  the  pressure  between  the  momentum  equation 
and  the  energy  equation  arises,  regardless  of  the  order  (in  u)  of  the  equilibrium  distribution  function  —  this 
is  true  even  in  the  continuum  case.  Furthermore,  certain  terms  in  the  pressure  tensor,  P^-,  were  arbitrarily 
omitted  in  the  Navier-Stokes  equation  derived  from  the  model  [62,  61].  Therefore,  the  use  of  P ij  is  ad  hoc 
at  best.  Following  the  above  analysis,  one  can  also  conclude  that  the  multi-component  model  constructed 
in  the  same  manner  [51,  50]  suffers  from  the  same  problems.  One  distinctive  feature  of  this  model  is  that, 
by  using  Cahn-Hilliard’s  model,  the  surface  tension  can  be  controlled  independently  of  the  equation  of  state 
by  the  density  gradient.  This  appears  to  be  the  reason  why  the  spurious  mass  flux  is  reduced  in  this  model 
[62,  61]. 

Most  likely,  the  aforementioned  problems  were  not  discovered  in  [62,  61,  51,  50],  because  all  the  tests 
conducted  in  [62,  61,  51,  50]  were  not  designed  to  test  hydrodynamics  per  se ,  they  were  designed  to  test  other 
properties  such  as  the  equilibrium  density  profiles  and  Laplace’s  law,  which  are  related  more  to  spinodal 
decomposition  than  to  hydrodynamics,  and  can  be  well  produced  by  other  lattice  Boltzmann  models  [21,  18, 
19,  20,  46],  or  even  models  without  hydrodynamics,  such  as  Cahn-Hilliard’s  model  [3,  13].  In  fact,  it  has  been 
shown  that  the  model  of  Swift  et  al.  does  not  satisfy  the  Navier-Stokes  equation,  in  either  the  bulk  region 
or  interface  region,  although  it  does  possess  conserved  quantities  [31].  Numerous  simple  hydrodynamic 
tests  showed  that  the  departure  of  the  model  from  the  Navier-Stokes  equation  is  quite  significant  both 
quantitatively  and  qualitatively  [31].  For  instance,  when  the  model  is  used  to  simulate  a  simple  hydrodynamic 
problem  such  as  a  single  droplet  subject  to  shear  [64],  the  simulation  [64]  does  not  produce  any  credible 
quantitative  hydrodynamic  results  because  the  model  does  not  even  satisfy  the  Navier-Stokes  equation.  The 
best  that  the  kind  of  simulations  presented  in  a  recent  work  by  Wagner  and  Yeomans  [64]  can  provide  is 
some  vague,  poor,  qualitative  description  of  flow  phenomena,  which  is  often  incorrect.  Therefore  we  doubt 
the  validity  of  numerical  predictions  made  by  the  model. 

It  should  also  be  noted  that  the  Hamiltonian  approach  [59,  60]  and  the  free  energy  approach  [62,  61,  51, 
50]  are  indeed  equivalent  in  terms  of  phenomenology.  Given  a  Hamiltonian  7~L  of  an  interacting  7V-particle 
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system,  the  corresponding  free  energy  ^  can  be  obtained  via  the  partition  function  based  upon  %: 


(5.11) 


N 


n  =  H 

i—  1 


+  U(Xi) 


i<j 


where  &  and  x\  are  the  phase  space  of  the  i-the  particle,  nii  is  the  particle  mass,  U{xi)  is  an  external  field 
and  V(\xi  —  Xj |)  is  a  mean-field  two-body  interaction  potential,  the  partition  function  is 

(5.12)  Z  =  J  dxd£  exp(— l-L/ksT) , 

where  (x,  £)  represents  the  entire  phase  space  of  the  7V-particle  system.  Consequently,  the  free  energy  is 
given  by 


(5.13) 


$  =  -kBT\nZ. 


Thus,  in  principle  information  is  neither  gained  nor  lost  whether  the  problem  is  formulated  in  terms  of  % 
or  One  cannot  claim  that  using  the  free  energy  and  utilizing  the  Maxwell  construction  leads  to  a  better 
or  more  physical  model.  Perhaps  the  only  advantage  of  using  the  free  energy  is  that  ^  is  a  global  state 
variable  and  therefore  it  is  independent  of  coordinates.  However,  this  advantage  bears  no  relevance  in  the 
LBE  models.  In  light  of  the  thorough  analysis  provided  above,  we  conclude  that  the  claim  made  by  Swift 
[62,  61]  et  al.  that  “we  show  for  the  first  time  that  it  is  possible  to  set  up  a  lattice  Boltzmann  scheme  modeling 
isothermal  hydrodynamics  for  two-phase  systems ”  is  unfounded. 

In  summary,  the  main  difference  between  the  model  derived  from  the  Enskog  equation  and  the  existing 
ones  is  in  the  physics.  Our  starting  point  is  the  Enskog  equation  for  dense  gases  in  which  the  non-ideal  gas 
effects  are  naturally  considered,  whereas  in  all  other  existing  models  [59,  60,  62,  61,  51,  50,  27],  the  starting 
point  is  the  original  Boltzmann  equation  which  is  only  suitable  for  dilute  gases  (ideal  gases) .  Beginning  with 
incorrect  physics,  one  has  no  choice  but  to  invoke  various  ad  hoc  approximations.  Inevitably,  various  defects 
exist  in  the  aforementioned  models.  One  notable  feature  common  to  these  models  is  that  the  viscosity  is 
independent  of  non-ideal  gas  effects,  which  is  inconsistent  with  the  Enskog  equation. 

6.  Conclusion.  We  are  now  in  the  position  to  lay  out  the  procedure  for  constructing  a  thermody¬ 
namically  consistent  lattice  Boltzmann  equation  for  non-ideal  gases.  Given  inter-particle  interactions,  the 
radius  distribution  function  g  can  be  computed  in  principle,  and  the  collision  term  in  the  Enskog  equation, 
J',  can  be  constructed.  This  collisional  term  would  correctly  produce  the  non-ideal  gas  effects.  With  this 
term  implemented  in  the  lattice  Boltzmann  model,  a  thermodynamic  and  hydrodynamic  consistency  can  be 
achieved  in  the  sense  of  the  finite  difference  approximation.  With  g  given,  the  free  energy  density,  can  be 
obtained  explicitly.  Subsequently,  other  pertinent  thermodynamic  quantities  such  as  the  equation  of  state, 
pressure  tensor,  surface  tension,  and  so  on,  can  be  directly  and  easily  derived  from  the  free  energy  while  the 
correct  hydrodynamics  is  preserved  in  the  lattice  Boltzmann  equation. 

It  should  be  emphasized  that  because  the  Boltzmann  equation  describes  mesoscopic  dynamics,  the 
constraints  imposed  on  it  must  be  compatible  with  the  mesoscopic  dynamics.  Specifically  speaking,  given 
an  arbitrary  interaction  in  a  system,  one  can  in  principle  compute  the  equation  of  state  (e.#.,  by  means  of 
the  virial  expansion).  This  is  an  averaging  process,  because  macroscopic  observables  (the  equation  of  state, 
surface  tension,  etc.)  are  averaged  macroscopic  quantities.  But,  the  reverse  are  not  true  in  general:  given 
an  arbitrary  equation  of  state,  one  may  not  be  to  able  to  uniquely  find  a  corresponding  interaction  in  the 
mesoscopic  description.  However,  this  can  be  achieved  in  the  formalism  of  the  lattice  Boltzmann  equation, 


16 


owing  to  the  simple  structure  of  the  formalism.  Nevertheless,  the  simplicity  of  the  LBE  method  does  not 
comes  without  a  price  —  the  inconsistency  in  the  LBE  thermodynamics  and  discrete  effects  are  inherent  to 
the  LBE  models.  In  contrast  to  the  LBE  method,  it  is  also  worth  while  to  mention  some  new  approaches 
more  closely  related  to  standard  computational  fluid  dynamics  methods  [34,  49,  33]  which  show  promise  for 
dealing  with  the  interfacial  problems. 

It  is  a  fair  observation  that  so  far  a  large  part  of  the  lattice  Boltzmann  enterprise  rests  upon  the  phe¬ 
nomenology  of  creating  various  “new”  equilibrium  distribution  functions  to  accommodate  different  physical 
phenomena  ranging  from  non-ideal  gases  or  multi-component  fluids  to  visco-elastic  media  [55].  Previous 
procedures  to  construct  the  equilibrium  distribution  can  be  summarized  as  follows.  By  observing  the  hy¬ 
drodynamic  equations  of  a  system  of  interest,  one  can  anticipate  those  terms  in  the  equilibrium  distribution 
which  are  necessary  to  produce  the  desired  results  (usually  a  desirable  stress  tensor).  Then  proportionality 
factors  for  these  terms  are  determined  by  the  conservation  constraints.  It  is  evident  that  this  approach  lacks 
mathematical  rigor  and  that  the  models  derived  in  this  fashion  may  suffer  from  artificial  defects  which  are 
uncontrollable,  such  as  the  models  in  Refs.  [62,  61,  51,  50,  55].  The  problem  common  to  these  models  is 
that  the  mathematical  rigor  of  the  Chapman-Enskog  analysis  has  been  completely  ignored,  as  typified  by 
the  work  in  [62,  61]. 

It  is  important  to  point  out  that  the  rigor  of  the  Chapman-Enskog  analysis  can  be  retained  without 
following  the  viewpoint  of  deriving  lattice  Boltzmann  models  via  discretization  of  the  corresponding  contin¬ 
uous  kinetic  equation.  Given  a  set  of  discrete  velocities  on  a  lattice  space  with  a  collision  operator  obeying 
conservation  laws  and  associated  symmetries,  an  orthogonal  basis  spanned  by  the  eigenvectors  of  the  collision 
operator  can  be  obtained  [8,  16,  17].  The  kinetic  modes  of  the  basis,  which  are  fluxes,  can  have  different 
relaxation  times  [8,  16,  17].  Not  only  does  this  approach  overcome  some  shortcomings  of  the  single  relaxation 
time  method  such  as  a  fixed  Prandtl  number,  but  it  also  follows  the  Chapman-Enskog  analysis  rigorously. 
We  suspect  that  this  approach  may  be  equivalent  to  a  discrete  version  of  the  hierarchy  of  kinetic  equations 
proposed  by  Levermore  [39]. 

In  this  paper  we  carry  out  a  systematic  derivation  of  the  lattice  Boltzmann  equation  for  non-ideal  gases 
from  the  Enskog  equation.  It  should  be  stressed  that  the  procedure  illustrated  here  is  general  and  can  be 
easily  extended  to  other  lattice  Boltzmann  models,  e.g .,  multi-component  models  [21,  18,  19,  20].  This 
procedure  can  also  be  used  to  improve  the  accuracy  of  the  lattice  Boltzmann  models  systematically.  Our 
procedure  can  be  briefly  summarized  as  follows.  First  of  all,  one  can  observe  the  equation  of  state  of  a  system, 
and  extract  the  non-ideal  gas  part  in  it.  This  part  is  related  to  the  radial  distribution  function  g.  From  g , 
the  additional  collision  term  responsible  for  it  can  be  constructed.  Then  one  can  systematically  discretize  the 
Enskog  equation  with  the  additional  collision  term  to  obtain  the  corresponding  lattice  Boltzmann  equation. 
This  approach  is  not  only  rigorous,  but  also  systematic.  The  equilibrium  distribution  is  uniquely  determined 
in  this  procedure.  It  enables  one  to  see  clearly  what  approximation  is  made  in  the  derivation  of  the  lattice 
Boltzmann  equation.  In  this  way  it  can  be  shown  that  the  accuracy  of  the  lattice  Boltzmann  equation  is 
indeed  first  order  in  time  and  second  order  in  space  [67]. 

In  additional  to  an  a  priori  derivation  of  the  lattice  Boltzmann  model  for  non-ideal  gases,  we  explicitly 
illustrate  the  differences  between  our  model  and  existing  ones.  Based  upon  our  analysis,  we  can  conclude 
that  the  problem  in  the  model  of  interacting  potential  [59,  60]  is  a  minor  one  and  can  be  easily  fixed  by 
either  directly  using  a  forcing  term  as  we  proposed,  or  adding  a  correction  to  remove  the  dependent 
terms.  In  contrast  to  the  model  of  interacting  potential,  the  model  of  free  energy  [62,  61]  presents  many 
major  problems.  It  starts  with  an  intention  to  correct  the  thermodynamic  inconsistency  in  other  models, 
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but  it  ends  up  with  more  serious  inconsistencies,  because  the  pressure  tensor  in  the  model  is  constructed 
without  any  physical  basis  at  the  level  of  the  Boltzmann  equation.  It  is  also  important  to  point  out  that 
none  of  these  models  can  lead  to  the  correct  energy  balance  equation,  and  therefore  they  are  inconsistent 
with  their  continuous  counterpart  —  the  Boltzmann  equation.  Starting  with  the  Enskog  equation  in  the 
presence  of  an  external  field  and  through  a  rigorous  discretization  procedure,  we  can  obtain  a  consistent 
thermodynamics  and  hydrodynamics  for  non-ideal  gases  in  the  sense  of  the  discretizing  approximation. 
With  this  systematic  means,  one  can  use  either  an  interaction  or  a  free  energy  to  obtain  the  equation  of 
state,  that,  when  incorporated  into  a  collisional  term,  accounts  for  non-ideal  gas  effects  among  the  particles. 

Our  future  work  will  extend  our  theory  to  multi-component  fluids,  and  obtain  a  consistent  thermody¬ 
namics  for  lattice  Boltzmann  models. 
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Appendix  A.  Modified  Boltzmann  Equation  for  Dense  Gases. 

The  Boltzmann  equation, 

(A.l)  dtf  +  |-V/  +  a-V^f  =  J , 

can  be  modified  for  dense  gas  by  explicitly  considering  the  volume  exclusion  effect  in  the  collision  term  J  for 
hard  spheres  of  radius  ro  as  the  following  [5,  30,  23]: 

(A. 2a)  J  =  J  d/j, i  [g(x  +  r0r)f(x,  £')/(*  +  2ro  f,  £|)  -  g(x  -  r0r)f(x,  £)f(x  -  2 r0r,  &)] 

=  J(0)  +  J(1)  +  J(2)  +  •  •  •  , 

(A. 2b)  J(0)  =g  Jdm(ffi~ffi), 

( A.2c)  J^=r0J  dm  (ff[  +ffi)f-Vg, 

(A. 2d)  J<2>  =  2 r0g  j  dm r  •  (/'V/J  +  /V/i) , 

where  g  is  the  radial  distribution  function,  r  is  the  unit  vector  pointing  from  the  particle  /i  to  the  particle 
/,  j(n\  n  —  0,  1,2,  ...,  are  obtained  by  the  Taylor  expansion  of  g ,  /i  =  /(|  i),  and  f[  =  f(£[)  in  Eq.  (A.  2a) 
of  J  about  x, 

(A. 3)  dfi\  =  |||i  -  |||crdnd|i  , 

and  a  and  ft  are  differential  collision  cross  section  and  the  solid  angle  in  coordinate  x-space.  The  Enskog 
equation  is  also  called  the  modified  Boltzmann  equation  for  dense  gases  in  the  literature  [23]. 

The  term  J given  by  Eq.  (A. 2b)  is  the  usual  collisional  term  in  the  Boltzmann  equation  with  an  extra 
factor  g ,  which  can  be  approximated  by  the  BGK  approximation,  i.e ., 

(A.4)  J( °)  =  -|[/ -/(")], 

where  /(0)  is  the  equilibrium  distribution  function  of  Maxwell  and  Boltzmann. 
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ro, 


The  term  and  J ^  can  be  explicitly  evaluated  for  hard  sphere  potential.  For  hard  spheres  of  radius 
we  have 


(A.5) 


(7  dCl  = 


ds  =  2ro  cos  tid'd  , 
sds  =  2 7q  sin(2$)d$d<p  , 


in  2D, 
in  3D, 


where  ti  is  the  azimuthal  angle  between  the  2- axis  parallel  to  (£1  —  $,)  and  (£'  —  £),  and  0  <  ti  <  tt/2,  tp  is 
the  polar  angle  on  the  plane  perpendicular  to  the  z-axis,  and 


(A.6) 


s  =  2ro  sin  ti 


is  the  impact  parameter  in  the  collision.  With  the  approximation  of  /  «  /(0)  in  Eqs.  (A. 2c)  and  (A. 2d),  we 
have 


(A. 7a) 


(A. 7b) 


jM  =  -f0)bp(t-u)-Vg, 


J{2)  =  -fi0)bpg 

+  {cdT2) 


2(£  —  u)-V\np  + 
'  (£-«)2 


(&  -  -  Uj)diUj 


0 


-  1 


(D  +  2) 

}v-u+i{ 


D 


0 

(d-uf 


(23  +  2)  0 


-1  (€ 


i/)-Vln0 


where  b  is  the  second  virial  coefficient,  p,  u ,  and  0  =  feeT/ra,  are  mass  density,  velocity,  and  normalized 
temperature  of  fluid,  respectively;  m  is  the  particle  mass,  and  D  is  the  dimension  of  the  £-space.  In  above 
equation,  the  Einstein  notation  for  summation  among  repeated  indices  is  used.  The  second  virial  coefficient 
for  the  hard  sphere  gas  is: 


(A. 8)  b  =  Vo /m  , 

where  Vo  is  the  volume  of  a  hard  sphere,  which  is  167rro/3  in  3D,  or  7 tTq  in  2D.  Note  that  bp  =  Vo n,  where 
n  is  the  particle  number  density,  is  a  dimensionless  quantity  called  packing  fraction  (for  the  hard-sphere 
system),  and  g  is  a  function  of  bp. 

It  should  be  noted  that  the  collision  term,  J,  in  the  Enskog  equation  does  not  conserve  mass,  momentum, 
and  energy  locally ,  because  it  involves  non-local  interactions  [5].  However,  mass,  momentum,  and  energy 
are  conserved  globally.  The  non-local  interaction  is  expected  to  produce  non-ideal  gas  effects  due  to  the 
exclusive  volume  in  momentum  and  energy  equations  in  hydrodynamics.  The  first  term  in  the  expansion  of 
J,  J(°)  ,  is  the  usual  collision  term  in  the  original  Boltzmann  equation  for  dilute  gases  (multiplied  by  a  factor 
g),  and  it  conserves  mass,  momentum,  and  energy  locally.  Other  terms,  J ^  for  n  >  1,  do  not  conserve 
mass,  momentum,  and  energy  locally,  they  are  responsible  for  the  flux  (of  mass,  momentum,  and  energy) 
transfer  due  to  non-local  interactions. 

To  the  first  order  approximation  in  the  Chapman-Enskog  analysis,  only  j(°),  J^l\  and  J ^  shall  be 
retained  in  the  modified  Boltzmann  equation  for  the  dense  gases.  Higher  order  collisional  terms,  j(n\  n  >  3, 
are  neglected  because  they  are  involved  higher  order  or  higher  power  of  gradients  of  p,  u ,  and  0.  The  term 
j(2)  can  be  simplify  for  incompressible  and  isothermal  fluids.  In  that  case,  the  last  two  terms  in  Eq.  (A. 7b) 
vanish.  Then  the  term  of  diUj  must  be  neglected  owing  to  the  conservation  constraints.  The  term  of  Vp  can 
be  included  into  J'  by 

(A. 9)  J'  =  -/(0)  bp  (£-«)•  [Vg  +  gV  In p2]  -  /<°>  bpg{Z~  u)-V\n{p2g) . 

It  is  clear  that  J'  conserves  mass  locally.  However,  it  is  responsible  for  flux  transfer  due  to  the  non-local 
interaction.  The  first  order  moment  of  J',  which  is  the  mass  flux,  gives  the  part  of  the  equation  of  state 
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attributed  to  non-ideal  gas  effects: 


(A. 10)  Jd££Jf  =  -bpg  J /(0)£  (£  —  u)  -  V  \n(p2 g)  -  b0p2gV\n(p2g)  = -V(0bp2g)  . 

Combining  with  the  ideal  gas  part  of  the  equation  of  state,  pO ,  we  obtain  the  total  equation  of  state: 
(A. 11)  P  =  p0(l  +  bpg). 

For  hard  sphere  gases,  the  radius  distribution  function  g  is  known  as  up  to  (bp)3  [5,  30]. 

The  energy  transfer  due  to  the  effect  of  J'  is: 


(A.12) 


\  fd££2J'  =  —^bpg  [d€f0)£,2(€-u)-X7ln(p2g)  -  u-X7(0bp2g) . 


This  correctly  corresponds  to  the  non-ideal  gas  thermal  equation  of  state,  Eq.  (A.  11). 

Appendix  B.  Chapman-Enskog  Analysis  of  Lattice  Boltzmann  Equation. 

By  introducing  the  following  expansions  [41,  24]: 


(B.la) 

(B.lb) 

(B.lc) 


fa(x  +  ea5t,  t  +  5t)  =  t). 

n= 0 

OO 

n= 0 
oo 


where  e  =  St  and  Vt  =  [dt  +  ea-V),  we  can  rewrite  the  lattice  Boltzmann  equation  with  a  forcing  term 
(B.2)  fa(x  +  eaSt,  t  +  8t)~  fa(x,  t)  =  -%«(*,  t)  -  f^(x,  t)}  +  J'Jt  -  Fa8t 

T 

in  the  consecutive  order  of  the  parameter  e  as  follows: 

(B.3a)  e°  :  , 

(B.3b)  e1  :  /«  =  -- Vtof  (°>  , 


-0  . 

m  - 

f  (eq) 

J  a 

J  a  ? 

.1  . 

fP  = 

--vt 

9 

2  . 

fi2)  = 

T  , 

J  OL 

“2 g  1 

where  Vtn  =  (dtn  +  ea-V).  Note  that  both  Fa  and  J'a  in  Eq.  (B.2)  do  not  appear  in  the  above  equations. 
However,  they  will  appear  in  the  governing  equations  in  what  follows.  The  distribution  function  fa  is  the 
normal  solution  which  is  constrained  by: 


(B.4a)  £/<?>=", 

a  e«  J  Pu  _ 

(B.4b)  £/<”)  1  =  °  n>  0, 

a  e«  J  0  J 

where  the  equilibrium  / ^  (for  the  9-bit  model)  is  given  by: 

CD  ^  Aeq)  J\  ,  o(e«M)  ,  9(ea-u)2  3  U2 

(B-5)  fa  q;  =«>*#>  0  +  3-^—  +  -— ^ - Y& 
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The  additional  collision  term,  J'a,  which  is  responsible  for  the  volume  exclusion  effect  in  dense  gases,  is  given 
by 


(B.6)  J'a  =  b  pg  (ea  —  u)  -V  \n(p2g) . 

The  forcing  term,  Fa,  for  the  nine- velocity  model,  is  given  by, 

1  (s  • u ) 

(B.7)  Fa  =  —2>wap  —r (ea  -u)  +  3ya  J  ea  a, 

cl  c4 

which  satisfies  the  following  constraints: 

(B.8a)  EF“  =  0’ 

a 

(B.8b)  y^eaFa  =  -pa, 

a 

(B.8c)  E  eajieajFa  =  -p  ( UiOj  +  ujch)  . 

a 

Also,  fa  is  a  Chapman-Enskog  ansatz ,  f(x,  t)  =  f(x,  p,  u ,  0),  i.e .,  the  temporal  dependence  of  fa  is 
through  the  hydrodynamic  variables  p,  u  and  0.  Therefore, 


(B-9) 

for  isothermal  fluids. 

For  the  9-bit  model,  \ 

^e  have: 

dtfa  =  ^dtP  +  ^dtu 

(B.lOa) 

E  1  =  p  > 

Q/ 

(B.lOb) 

Ee“i«0)  = 

Cf 

(B.lOc) 

E  ea,iea,jfa)  =  , 

Cf 

(B.lOd) 

a 

and 

(B.lla) 

EJ«  = 

OL 

0, 

(B.llb) 

eaJa 

oc 

=  -9bX7(p2g), 

(B.llc) 

Ee«he' 

a,jj'a  =  b  [ UiUjU-V  -  6  (Uidj  +  Ujdi )]  (p2p) 

a 


where  Sij  and  Sijki  are  the  Kronecker  delta  with  two  and  four  indices,  respectively,  and 

(B.12)  A ijkl  —  dijdjd  T  dik^ji  T  du^jk  • 


The  governing  equations  of  fa  up  to  the  order  of  e  are: 

(B.13a)  ©t0/(0)  =  -£/(1)+^--Fa, 

T 

(B.13b)  +  (2T~  g)A0/(1)  =  -£/(2)  - 

ZT  t 
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In  the  derivation  of  Eq.  (B.13b),  we  have  made  the  following  approximation  that  Vtog  ~  0,  which  is  accurate 
up  to  0{u2).  The  moments  of  the  first  order  governing  equation,  Eq.  (B.13a),  lead  to  the  Euler  equations: 

(B.14a)  dtop  +  V-(pu)  =  0, 

(B.14b)  dto(pu)  +  V-n(°)  =  F  -0bV(p2g) , 

where  n<°)  =  J2aec*eafa)  is  the  zeroth-order  momentum  flux  tensor.  With  11 given  by  Eq.  (B.lOc),  the 
above  equations  can  be  rewritten  as: 


(B.15a)  dt0P  +  V-(pu)  =  0, 

(B.15b)  dtou  +  u  Vu  =  — -VP  +  a  , 

where  a  =  F / p  is  the  acceleration,  and 


(B.16)  P  =  p0(l  +  bpg) 

is  the  equation  of  state  for  non-ideal  gas,  depending  on  the  radial  distribution  function  g.  (Note  that  for  the 
nine- velocity  isothermal  model  here,  0  =  c2/ 3.) 

The  moments  of  the  second  order  governing  equation,  Eq.  (B.13b),  lead  to  the  following  equations: 


(B.17a)  dtlp  =  0, 

(B.17b)  dtl(pu)+  (2r~g)V-nW  =0, 

At 

where  nW  =  J2a  eaeafa  ^  is  the  first-order  momentum  flux  tensor.  With  the  aid  of  Eqs.  (B.10)  and  (B.15), 
we  have: 


ni)  =  =  — ^2eaiieajVto  /^0)  =  —  [dt0^ij  +0(V-pu5ij  +  dipuj  +  djpUi) 

a  9  a  9 

=  -T-  [0  (dtop  +  V-pu)  Sij  +  dt0  (puiUj)  +  0  ( dipuj  +  djpui)]  =  -^0  p(diUj  +  djUi)  +  0(M3) , 


where  di  =  d/dxi.  In  the  above  result  of  11 the  terms  such  as  Uidjp  have  been  neglected  because  Vp  is 
of  the  order  0(M2),  and  it  is  done  consistent  with  the  small  velocity  expansion  of  fa  up  to  the  order  of 
0(u2).  [Note  that  0(u)  =  0(M ),  therefore  we  take  the  liberty  to  interchange  these  notations.]  Similarly, 
we  have 


djU^  =  -^<9  pdj(diUj  +djiii)  +  0(M3) 

(B.18)  =  -^6 p(diV-u  +  V2Ui)  +  0(M3)  =  -T- 6pV2Ui  +  0(M 3) , 

where  the  term  V"U  has  been  neglected  because  it  is  of  0(M2)  due  to  Eq.  (B.15a). 

Combining  the  first  and  the  second  order  results  [Eqs.  (B.14)  and  (B.  17)]  together  by  dt  =  dt0  and 

recalling  that  e  =  St,  we  have  the  Navier-Stokes  equations  [accurate  up  to  the  order  of  0(M2)  in  momentum 
equation] : 

(B.19a) 

(B.19b) 
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dtp  +  V-(pu)  =  0, 

dtu  +  u-Vu  =  — -VP  +  vSJ2u  +  a  , 
P 


where  the  viscosity  is  given  by 


„  =  (L-1)))St=(lLZ&f, 

9  2  6 g  St 


and  the  pressure  (the  equation  of  state)  is  given  by 
(B.20) 

where  0  =  c2/ 3  has  been  substituted.  With  the  above  equation  of  state,  the  sound  speed,  cs,  is  given  by 


P  =  p8{l  +  bpg)  =  -c2  p{l  +  bpg), 


(B.21) 


c2  =  9 


3° 


1 + 


It  should  be  pointed  out  that,  if  instead  of  Eq.  (B.8c),  the  following  constraint  is  imposed: 


(B.22)  ^  ea,ieaj-Fa  =  0 , 

a 

then,  the  term  pa  u ,  which  is  the  work  done  by  the  force,  does  not  appear  in  the  energy  balance  equation. 
Therefore,  the  constraint  of  Eq.  (B.8c)  must  be  imposed  to  assure  a  correct  energy  balance  equation. 

Appendix  C.  Equilibrium  Distribution  Function  Shifted  by  Acceleration. 

If  we  start  with  the  BGK  Boltzmann  equation  without  a  forcing  term: 

(c.i)  dtf+^f  =  -\\f-r o)], 

T0t 

and  assume  that  particle  is  impulsively  accelerated  by  acceleration  a  with  the  mean  free  time  rSt-  Under 
this  circumstance,  the  equilibrium  distribution  function  becomes  [66]: 


(C.2) 

Accordingly, 


f{0)(p,  u  -  arSt,  9)  =  p  (2w9)  D /2  exp  [-(£  -  u  +  ar<5t)2/2i9]  . 


(C.3) 
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Here  we  have  consistently  ignored  the  terms  of  the  order  0(S2)  or  higher  order.  Substituting  0  =  c2/ 3,  where 
c  =  Sx/ St,  we  have 


(C.4)  f(^=waP 


1  + 


3(ea  ■  u)  9(ea  ■  u)2  3 u2 
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2c4 


2c2 


-3  wap  rSt 


1  3 (ea-u) 

~o  \e&  ~  U)  4 - 4 -  eo 


• a . 


The  second  part  of  fa exactly  produces  the  forcing  term  Fa  obtained  previously  when  fa  is  substituted 
into  the  following  lattice  Boltzmann  equation  without  a  forcing  term: 


(C.5) 


fa(x  +  ea5t,  t  +  5t)  —  fa(x,  t)  =  --  \fa(x,  t)  -  f(aeq)(x,  t ) 

T  L 
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